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Abstract. We report a new analytical method for exact solution of homoge- 
neous linear ordinary differential equations with arbitrary order and variable 
coefficients. The method is based on the definition of jump transfer matrices 
and their extension into limiting differential form. The approach reduces the 
nth-order differential equation to a system of n linear differential equations 
with unity order. The full analytical solution is then found by the pertur- 
bation technique. The important feature of the presented method is that it 
deals with the evolution of independent solutions, rather than its derivatives. 
We prove the validity of method by direct substitution of the solution in the 
original differential equation. We discuss the general properties of differential 
transfer matrices and present several analytical examples, showing the appli- 
cability of the method. We show that the Abel-Liouvillc-Ostogradski theorem 
can be easily recovered through this approach. 



1. Introduction 

Solution of ordinary differential equations (ODEs) is in general possible by dif- 
ferent methods [1]. Lie group analysis [2-6] can be effectively used in cases when 
a symmetry algebra exists. Factorization methods are reported for reduction of 
ODEs into linear autonomous forms [7,8] with constant coefficients, which can be 
readily solved. But these methods require some symbolic transformations and are 
complicated for general higher order ODEs, and thus not suitable for many com- 
putational purposes where the form of ODE is not available in a simple closed 
form. 

General second-order linear ODEs can be solved in terms of special function [9] 
and analytic continuation [10]. Spectral methods [11,12] and orthogonal polynomial 
[13] solutions of linear ODEs have been also reported. Solution properties of linear 
ODEs with arbitrary orders are discussed in a report [14] , and several approximation 
methods have been also published [15]. All of these methods express the solution in 
terms of infinite series, relying mainly on numerical computations. Analytic solution 
of such equations has been also recently presented [16], but only at a singular point. 

Linear ODEs of the nth-order can also be transformed to a system of n linear 
first-order differential equations by reduction method [17], where the derivative of 
one parameter is regarded as the next parameter. This approach is based on the 
evolution of the solution function and its derivates up to the order n. This method 
is therefore not useful when the evolution of independent solutions is of interest 
rather than the derivatives, such as waves in inhomogeneous media. 
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Recently, we have developed a new matrix-based exact analytical method in 
order to study the propagation of electromagnetic waves in non-uniform media 
with an arbitrary refractive index function [18,19]. The method is based on the 
modification of the well-known transfer matrix theory in optics [20] and quantum 
mechanics [21], and its extension into differential form. This method is simple, 
exact, and efficient, while reflecting the basic properties of the physical system 
described by the differential equation. We have also developed proper formalism 
for anisotropic non-uniform media, in which the governing equations are no longer 
scalar [22]. In this method, the solution is found by integration of an exponent 
square matrix and then taking its matrix exponential. 

In this paper, we show that the differential transfer matrix method (DTMM) 
[18,19] can be also used for exact solution of linear ODEs with arbitrary variable 
coefficients. We only consider the homogeneous equation and its linear indepen- 
dent solutions, since once the linear independent solutions are known, the particular 
solution can be found by the method of variation of parameters due to Lagrange 
[2] . We show that the nth order ODE is transformed into a system of n homoge- 
neous first order linear equations, which can be integrated either numerically by 
Runge-Kutta [23,24] or reflexive [25] methods, or analytically by matrix exponen- 
tial techniques. Through direct substitution we rigorously show that the presented 
analytical solution satisfies the corresponding differential equation, and justify the 
approach through several illustrative examples. We also show that the famous 
Abel-Liouville-Ostogradsky theorem for the Wronskian of a linear ODE [26-28] can 
be easily proved through the DTMM. 

In §2 we describe the outline of the problem and then formulate the jump transfer 
matrix method in §3. In §4 differential formulation, the fundamental theorem of 
DTMM, and treatment of singularities are described. Finally, the application of 
the DTMM to several ODEs of various orders is presented. 



2. Statement of Problem 
Consider the homogeneous ordinary linear differential equation of order n 



(1) L/(s) = 0, /:S^C 

where / (x) is an unknown analytic complex function in the set of complex variables 
C with the connected domain § C C, and L is a linear operator given by 



(2) L = ^ a m (x) , a n (x) = l. 

m=0 

Here we assume that a m :§nC,m = 0, . . . ,n— 1 are arbitrary analytic complex 
functions of x, being referred to as variable coefficients. If coefficients are constant, 
then under the condition described just below, a general solution to JIJ is [1,29] 



(3) 



/ (x) = a exp {k\x) + c 2 exp (k 2 x) + 



+ c„ exp (k n x) , 
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where Cj, i = 1, . . . , n are complex constants determined from either boundary or 
initial conditions, and ki,i = l,...,n, being referred here to as wavenumbers, 
satisfy the characteristic equation 

n 

(4) ]Ta m fc™ = 0, i=l,...,n. 

Moreover, for © to be a general solution of with constant coefficients, it is 
necessary that Vi, j; h ^ kj, i.e. (@J has exactly n distinct roots. Now, if coefficient 
functions are not constant, then we propose a general solution to Q of the form 

(5) / (x) = fx (x) exp [hi (x) x] + f 2 (x) exp [k 2 (x) x] H h /„ (x) exp [k n (x) x] , 

or equivalently 

(6) f(x)=exp[&(x)} t F(x), 



' h (x) ' 




k\ (x) x 


h(x) 




k 2 (x) x 


, * (x) = 


_ fn (x) _ 




kn \X) X 



where the superscript t denotes the transposed matrix. Hereinafter, we shall refer 
to F (x) and $ (x) respectively as envelope and phase vector functions. In the 
above equation fi (x) are functions to be determined and ki(x), being referred to 
as wavenumber functions, are complex functions of x satisfying the generalized 
algebraic characteristic equation 



(8) J2 a - (*) ^ fr) = °' * = 1, 

For the sake of convenience, we here define the exponential of a vector v = 
{vi} nxl as exp (v) = {exp (vi)} nxl . In contrast, exp (M) when M is a square 
matrix represents the matrix exponential of M, given by 



(9) exp(M) =1+ — M ™ 

m— 1 

We define a domain § to be non-degenerate if © has exactly n distinct roots 
for all x £ S, and m-fold degenerate if it is non-degenerate everywhere, but at 
finite number of isolated points, being referred to as singularities, at which (JSJ) 
has a minimum of n — m distinct roots. If S is at least 1-fold degenerate for all 
x 6 8, then S is referred to as entirely degenerate. Later, we shall observe that non- 
singular solutions require § to be non-degenerate, and degenerate domains require 
special treatment. 

We also define the diagonal wavenumber matrix as 



(10) 



K (x) = [h (x)Sij] nXn = diag[fei(x), . . .,k n (x)} , 
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with 8ij being the Kronecker delta. Obviously, the diagonal elements of K (x), and 
thus K (x) itself must satisfy ©■ Moreover, |K(x)| = ao(x) and tr{K(x)} = 
a n _i (x), where tr {•} denotes the trace operator. The phase vector and wavenum- 
ber matrix are thus related as # (x) = xJi (x) l„ x i, where l nx i is a vector with 
unity elements. Hence, one has exp (x)] = exp [.iK (x)] l nx i. 

3. Jump Transfer Matrices 

3.1. Transfer Matrix of a Single Jump. Without loss of generality, suppose 
that §Cl with R being the set of real numbers. Furthermore, we let the variable 
coefficients di(x), i = 1, . . . , n to be such stepwise constant discontinuous functions 
at x — X that K (x < X) = Ka and K (x > X) = Kb become constant matrices. 
We furthermore accept that both the sub-spaces §a — {x\ x £ S, x < X} and §b = 
{x\ x £ §, x > X} are non-degenerate. In this case, the solution in § would be given 

by 



(11) 



j exp[$A (a:)]' Fa, x £ §a 
1 " ' 1 exp[* fl (a:)]* F Bj x £ § B 



where the envelope vectors Fa and Fb must be constant, and $a (x) = {Akix} nXl , 
d?B {x) = {Bkix\ ny l with Aki and ski satisfying the characteristic equation J2J in 
Sa and Sb- Analyticity of / (x) across x = X requires that 



(12) 



/M = (X+) , m = 0,...,n-l, 



where /( m ) (x) is the mth derivative of / (x) with respect to x. 
results in the set of equations 



Expanding 1(12)1 



(13) D B exp (IK B ) F b = B A exp (XK A ) F A , 

in which Ka and Kb are given in 1101) . Also Da — Ia^V 1 ] and D; 

' ' L J J nxn 

[b^ 1 ] are matrices of Vandermonde type [30] given by 



(14) 



D 



1 

Aki 



1 

Ak 2 



■~ n - 1 A k^- X 



i to- 

aK 

One can rewrite l|13|) as 



1 

Akr, 



1 



1 

Bk 2 



i n— 1 un—1 
B k 1 B K 2 



1 



(15) 



Fb — Qa^bF/ 



(16) 



Qa^b = exp (-IK B ) D^Da exp (XK A ) , 



in which Qa^b is defined as the jump transfer matrix from A to B. Similarly, one 
can define the downward jump transfer matrix as Qb^a = QaLb- 
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The determinant of the jump transfer matrix Qa^b is 

\Va\ 



(17) \Qa^b\ =exp[X(tr{K A }-tr{K s })] 



ID, 



in which we have used the identity |exp(M)| = exp(tr {M}). The determinants can 
be also evaluated by noting the fact that |D r | are determinants of Vandermonde 
matrices, given by [30] 

f 1, n = l 

11 \rh - r kj), n > 1 
Here, the product operator runs over all possible ordered pairs (i,j) with i > j. 

n 

Now since tr {K r } = ^ r ki = r a„_i with r = A,B, l|17|) can be expanded as 



(19) |Qa^b| = exp [X(Aa„_i - B a n -i)} 



(a^i — Akj) 

. (sfc* - Bkj) 



3.2. Transfer Matrix of Multiple Jumps. When variable coefficients a,i(x), i — 
1, ...,n are stepwise functions of x with arbitrary number of discontinuities or 
jumps over interfaces, K (x) can be expressed as 

(20) K(x) = K„, xeS p , p = 0,...,P, 

where the subset E> p = {x\ X p < x < X p+ i} ,p = 0, . . . ,P is referred to as the pth 
layer. Obviously, X — inf {§} = inf {§ } an d X P+ i = sup{§} = sup{Sp}. In this 
case, the corresponding envelope vectors are related as 

(21) F s = Q r ^ s F r , r,s = 0,...,P, 

where Q r ^ s is the transfer matrix from layer r to s, obtained by multiplication of 
single jump transfer matrices as 



(22) Qr— >s = Qs-1— >sQs-2->s-l - • • Qr- 



>r+l- 



3.3. Properties of Jump Transfer Matrix. The transfer matrix Q r ^ s satisfies 
the basic properties 

(23) Qr^r = I, (self — projection) 

(24) Q,^s = Q7-ir: {inversion) 

(25) Q r ^ s = Qu^sQr^u, (decomposition) 
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(26) |Q r 



exp 



_p=r i=l 



n 

J i>3 



(ski skj) 



(determinant) 



with r, s,u = 0, . . . , P. The properties Ij23(l . I|24|l . and l|25|l are direct results of the 
definition of transfer matrices, and the determinant property (|26|l follows (|22|l and 
P|l. 

From (|19|l it can be observed that for systems with vanishing r a„_i, r = 0, . . . , P, 
the determinant property simplifies into 



(27) l^^UTJ^y 

This expression has the important feature of being dependent only on the local 
properties of the initial layer r and target layer s. 

The transfer matrix Q r _> s satisfies the scaling property, that is it remains un- 
changed under transformations X p — > aX p and r ki — ► a r ki- This can be directly 
observed cither from Q by making the change of variable x — ► ax, which accord- 
ingly scales the roots in an inverse manner, or direct substitution in (|16|) . ft also 
preserves the shifting property 



(28) Qr^s = exp (-£K S ) Q r ^ s exp (£K r ) , (shifting) 

in which £ is the amount of shift over x-axis and Q r ^ s is the transfer matrix 
corresponding to shifted space. This property follows (|f 61) and (|22|l . The shift the- 
orem plays an important role in the theory of wave propagation in one-dimensional 
periodic structures in electromagnetic and optics [31]. 

In the next section we shall show that how one can solve Q in its most general 
form by extension of jump transfer matrices into differential form. 

4. Differential Transfer Matrix Method (DTMM) 

4.1. Formulation. We first assume that the wavenumber matrix K (x) is a smoothly 
varying function of x. For the moment, we also assume that S is non-degenerate. 
Following the notation in the previous section, let A and B denote respectively 
the sub-domains E>a = [x — Sx, x[ and Sb = [x, x + Sx[. If 5x is small enough, 
one can approximate the wavenumber matrix as K (t € §i) ~ K (x) = K^, and 
K (t £ Sb) ~ K (x + 5x) e Kg. Accordingly, one has 



(29) F (x + 8x) » F B = Qa^bFa « Qa^bF (x) 

The above equation is accurate to the first order. So, one can expand the transfer 
matrix Qa^s through its definition l|16|) as 



Qa^b = exp (~xK B ) B A exp (xK A ) 

= exp [-x (Ka + SK)] (Da + ST>)~ 1 T>a exp (xK A ) 
(30) w exp (-iK A ) (I - x(5K) [D^ 1 - D A 1 ((5D)D^ 1 ] D A exp (j:K a ) , 

= exp 1-xKa) (I - z<$K) (I - D^JD) exp (xKa) 
» exp (-sKa) (I - x6K - D^JD) exp (xKa) 
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in which we have used the Taylor expansion of matrix exponential and differential 
property of inverse matrices [32], and have neglected the 2nd- and higher-order 
variations. The matrix <5K can be approximated as 



(31) 



d~K = [SkiSij] « [kl (x) Sij] dx = — — K (x) 8x 



From (|14|l one also has 
(32) 



(5D = 



(n~l) A k?- 2 (n-l) A k 2 



= CSK 
So we can rewrite l|2*9"|) as 



dx 



■ • 
• • 1 



(n - 1) A k r r , 



£K(x)Sx 



(33) SF (x) m (Qa^b — I) F (x) = U (x) F (x) Sx. 

But the matrix U (x) can be found from l|3U|) and (|32|l as 



(34) 



§K 



U (x) « - exp (-xK) (xl + D _1 C) -— exp (xK) , 

ox 



where the trivial subscript A has been dropped. 

Now, if we let 5x approach zero, (|33J) and l|34(l become exact and can be rewritten 

as 



(35) 



dF (x) = U (x) F (x) dx, 



(36) U (x) = -xK' (x) - exp [-xK (x)] D (x)" 1 C (x) K' (x) exp [xK (x)] . 

Hereinafter, we refer to U (x) as the differential transfer or the kernel matrix. In 
(l3l)l) . we have 



(37) 



C (x) - [(< - 1) k)- 2 (x)l , D (x) = [fcj- 1 (x)l 



(38) 



K (x) = [h (x) %]„ x „ , K' (x) = ^ (x) %]„ x „ 



While l|35|l can be integrated by numerical methods [23-25], it permits an exact 
solution of the form 



(39) 



F(x 2 ) = Q Xl ^ 2 F(x!), 
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for Vxi,X2 G Here Q Xl —> X2 is referred to as the transfer matrix from xi to X2, 
given by the Dyson's perturbation theory as [33-35] 



(40) 



Q X1 ^ X2 = I + JdhUih) + Jdtt / dt 2 U(ii)U(i 2 ) - 

Xi Xi X\ 

Jdh j dt 2 J dt 3 ... / dt m V(h) . . . U(i m ) + 

X\ X\ X\ X\ 



By defining P as the so-called chronological ordering operator [33-35], it is possible 
to rewrite l|40(l as 



X2 X 2 X 2 



(41) Q X1 ^ X2 = J^^ J dhj dt 2 f dt 3 ... f dt m ¥ [U(ti) . . . 17(0] 

7n ~ ^ ^ xi xi 

Often, ||1T]> is symbolically written as [33-35] 



(42) 



Qa;i— *i 2 



"exp 



U (x) dec 



) exp(M ;cl ^ 2 ), 



in which M 2;i ^ K2 is the integral of the kernel matrix and referred to as the transfer 
exponent matrix. The above expression can be greatly simplified if the kernel matrix 
U(a;) satisfies one of the few existing sufficient conditions for integrability, including 
when it commutes with Mj.^^ [36-38], or it satisfies the Fedorov's condition [38], 
or it commutes with itself for Vxi,X2 € S. In this case, one has 



(43) 



Qrri 



exp 



U (x) dx 



Although in general (|43|) is not necessarily the exact transfer matrix, it is known 
that both H39fl and (|43|l must have the same determinant [39] . Also their traces are 
equal at least to the second order. This can be easily justified by comparing the 
expansion of l|43|l to 1)40(1 and using the cyclic permutation property of tr{-} [30]. 
This means for the case of second-order equations with n = 2, H43|l can be used 
instead of ((42JI . with an accuracy better than second-order. If needed, the exact 
transfer matrix can be directly found by numerical integration of equation [39] 



(44) dQ Xl ^ x = U (x) Q xi _ x dx, 

with the initial condition Q^^j,, = I. 

4.2. Properties of Transfer Matrix. The transfer matrix Q Xl ^ X2 from x\ to 
X2 clearly preserves the properties iffi^ . (|2"4")l and (jJEJ), and as discussed above, 
its determinant is always equal to the determinant of its counterpart given by 
0. Therefore, |Q 

xi— >2J2 1 — exp(tr ^M. xi _>g; 2 }). But the transfer exponent matrix 
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M I1 ^ l2 can be written as the sum of the jump J Xl —> X2 and propagation T Xl _» X2 
matrices given by 



(45) 



X2 

f d 
■ I exp [— iK (x)] D (x) -1 — {D (x) exp [xK (a;)]} dx, 
J dx 



(46) 



K (x) dx. 



The validity of the identity M ai _ 



>X2 < J X1^X2 



T X1 _> X2 can be verified by com- 



paring the integrands of both sides. The propagation matrix T X1 ^ X2 is diagonal 
and also has a diagonal matrix exponential given by 



(47) 



exp (T Xl -nc a ) = 



exp / ki (x) dx 5i 



Therefore, its determinant can be easily found by multiplying its diagonal elements. 
The determinant 1 3 Xl —t X2 | can be obtained using the relation D' (x) — C (x) K' (x) 
according to lj^2(l . and the identity (see Appendix A) 



(48) 



as 



exp 



■'2 

J H(a:)" 1 H' (x) dx 



= \n- 1 (x 1 )n(x 2 )\ 



(49) 



exp (J Xl ->x 3 )\ = exp [-x 2 K (x 2 )] D (x 2 ) 1 D (xi) exp [x%K (xi) 



Finally, the determinant (Q^^;^ | can be thus found using the well-known iden- 
tities [29] |exp(A + B)| = |exp(A)| x |exp(B)|, and |exp(M)| = exp (tr {M}) as 



(50) |Q | = exp (tr {T }) x exp^r-tJ^^J) 

Using the (|47|l . (|4*9"|) . and (fTHj) one can finally obtain 



(51) 

|Qa:i-^a:2 I 



= exp( a : 1 tr{K( aJl )}-x a tr{K( a : a )}-J ? tr{K( a ;)}d aJ ) x n fefepfefe 

\ x i / i > i 



exp 



x x a n -\ (xi) - x 2 a n -i (x 2 ) - J a„_i (a:) dx 



i>j 

x "Q fei(xi)-fe 3 (xi) 



This equation can be compared to (|26|l as the determinant property of transfer 
matrix. Notice that if a n _i (x) = 0, then (|51|l reduces to 



(52) 



|Qxx- 



>X2 I 



n 

i>j 



ki (xi) - fe 3 (xi) 
fej (x 2 ) - kj (x 2 )' 
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This shows that the determinant of the transfer matrix is only a function of starting 
and end points x\ and X2- Note that the transformation 



(53) / (x) —> exp [-i Ja n -i (x) dx] h (x) 

in Q always gives rise to a new nth-order ODE with identically vanishing coefficient 
of M™" 1 ) (x) term. Therefore, the above property of determinants can be always 
met through the abovementioned transformation. 

The shifting property (23) and scaling properties must also preserve, since the 
transfer matrix Q Xl ^ X2 is indeed obtained by dividing the sub-domain [x\, x%] into 
infinitely many thin layers and multiplying the infinitesimal jump transfer matrices. 

4.3. Justification of DTMM. Here we prove that the DTMM formulation is 
exact by showing that the solution satisfies Q through direct substitution. 

Lemma 1. Let the Vandermonde matrix D be invertible. Then the elements of 

n 

its inverse D _1 = [7^] satisfy ^2 lirk™ 1-1 = S mr where 1 < m < n. 

i=i 

Proof. The Vandermonde coefficients 7^ are found from the expansion of La- 
grange interpolation polynomials as [40,41] 



l,...,n. 

n 
i=l 

n n n / n \ n 

(55) kj- 1 = J2 kr 1 E = E E ^ = J2 

i— 1 r— 1 r— 1 \i=l / r— 1 

Given a fixed m, the above equation can be transformed into a linear system of 
equations as A*D = B*, with D being the Vandermonde matrix, A = [a rm ] being 
the column vector of unkowns, and B = being the input column vector. 

Since D is invertible by assumption, then A and thus a rm must be unique. However, 
(|55|l is satisfied through setting a rm = 5 rm for 1 < m < n, and hence the proof is 
complete. □ 

Lemma 2 (Derivative Lemma). Suppose that P (x) is an n x 1 vector func- 
tion, defined on the connected non-degenerate domain §, satisfying the differential 
equation l|35|) with the kernel matrix given by (|36|) . Then, we have 



(54) r.w^nfr^ 7 ^ 1 ' 1 

J 3=1 

n 

Obviously T, (kj) = £ ^Xf 1 = Therefore A:™ -1 



d m ( ~\ 
(56) — [exp [* [x)\ 1 F (x)j = exp [* (x)} 1 K m (x) F (x) , m = 0, . . . , n, 

where K (x) is the wavenumber matrix given by IjlUI) and 3? (x) is the corresponding 
phase vector defined in (JJJ . 

Proof. For m = the proof is trivial. Also for m>0we have 



^{expf*)'^} 

( 57 ) = £ exp (*) * K m F + exp (*) * £K m F + exp (#) * K m £F > 

= exp ($) * (K m+1 + xK m K' + mK m - 1 'K' + K m U) F 
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where the dependence of matrices on x is not shown for the sake of convenience. 
From the definition of the kernel matrix (|36J) . (|57|) can be simplified as 



£{exp (*)'K" 1 f} 



(58) 



exp (*) K 



m+li 



exp (*) 1 [mK" 1 -^' - K m exp (-iK) D^CK' exp (xK)] F » 
= exp (*)'K m+1 F + l^ xl (mK™- 1 - IC^D^C) K'exp(xK)F 
= exp(#) t K m+1 F + li x „VK'exp(xK)F 

where l nX i is defined under i|10|) and V = toK" 1-1 — K m D _1 C . Now we show 
that the row vector li x „V is identically zero, but only if m < n. By defining the 
elements of the inverse of the Vandermonde matrix as D 1 = pyy], we have 



lixnV = ml lxn K r ' 



(59) 



. 1Xn K m I>- l C 

n n 

E E7ir(r-l)feJ- 2 fe[' 

a—I r— 1 



771 k 



n—1 n 
r— 1 i—l 



Since the inverse of D must exist by assumption, we thus have according to Lemma 
1 



(60) 



llXnV = 



n-1 



= Oixnj m < n - 



Therefore, the proof of theorem follows immediately by induction. □ 

Theorem 1 (Fundamental theorem of differential transfer matrix method). 

Suppose that F (x) satisfies the derivative Lemma, with the kernel matrix K (x) sat- 
isfying the characteristic equation (JSJ. Then, the function f (x) defined by is a 
solution of the differential equation 

Proof. By expansion of the operator L given in J3J on / (x) given in ©, and 
using the derivative Lemma we have 



Lf (x) = £ a m (x) (x) 



m=0 



(61) 



E a m (x)^{exp[*(x)] 4 F(x)} 



m=0 

£ a m (x)exp[$( 2 ;)] t K m (x)F(x) 

m— 



exp [€> (x)] ' 



E a m (x)K m (a 



F(x) 



But the summation within the brackets is the null matrix because of © • Thus the 
right-hand-side of the above equation must be identically zero. This completes the 
proof. □ 

4.4. Linear Independent Solutions. Consider the set of vectors Fj £ §™,z = 
l,...,n, forming a basis in If Q is a non-singular matrix, then the vectors 
Gi = QF^ 6 C n ,i = 1, . .. ,n would clearly constitute another basis on § n . Let S 
be a non-degenerate domain with the differential equation and solution given by 
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(|39J) . Since § is non-degenerate, then Q Xl ^x 2 must be non-singular for all x\, xi £ 
§. Therefore the set of n vector functions defined by Gi (x) — Q a ^ x Fi (a) € §™ 
form a basis on Let \fr denote the set of all functions operating on §, as ^ = 
: 8 i-> C}. We define the scalar functions g t £ * as gi (x) = exp 

[* (*)] Gi (x) 

and show that they form a basis, by inspection of their Wronskian determinant. 

Theorem 2. With the above definitions, the set of functions gi (x) , i = 1, . . . , n 
form n linear independent solutions on ^ . 



gt l) (*) 



Using 



Proof. The Wronskian matrix in \P can be written as W = 

Theorem 1 one has W = [exp (3?) K J ~ 1 Q a ^. E Fi] , with the dependences on x omit- 
ted. Now we define the nxn matrix F = [Fj]. Then W = Dcxp(a;K) Q^jF. 
Because S is non-degenerate |D| ^ 0. Also, Q a _ ^ must be non-singular and thus 
IQa^zl 7^ 0. Now since according to the assumption, F^ form a basis on |F| ^ 
and therefore for the Wronskian determinant we have |W| =^ 0. Thus, the proof is 
complete. □ 

4.5. Treatment of Singularities. Now let S is degenerate at a finite number of 
isolated singularities € §,« = 1,...,£, at which |D(£.;)| = 0. Therefore, the 
U(£i) is singular and M. Xl ^ X2 does not exist. Without loss of generality, we can 
assume S = 1. If S represents the integration domain [xi,^], the total transfer 
matrix Q Xl ^ X2 by the decomposition property (|25|) can be written as 



(62) Q^i — >X2 Q^+<5.T — txzQ^ — Sx — >£-\-6x Q:ei — >£ — Sx- 

The transfer matrices Qa^—^-sa, and Q^-\-s x ^x 2 do exist and can be evaluated di- 
rectly by integration and exponentiation of U(ir), as long as Sx is positive and 
finite. The transfer matrix Q£-5x->£+6x enclosing the singularity, also can be ap- 
proximated by its equivalent jump transfer matrix as 

/f>n\ Qn-8x->t;+6x ~ exp [- (£ + Sx) K (£ + Sx)] D 1 (£ + Sx) x 

( > D (£ - Sx) exp [(£ - Sx) K (£ - Sx)] 

This approach permits evaluation of the total transfer matrix Q Xl ^ X2 by making 
finite jumps over singularites. 

If § is entirely degenerate, then the transformation / (x) — > w (x) h (x) with 
to (x) being a non-constant function, results in a completely different characteristic 
equation (jSJ. A proper choice of w(x) leads to a new §, which can be no longer 
entirely degenerate. Then DTMM can be used to solve for h (x). 

5. Examples 

5.1. First-order ODEs. It is easy to show the consistency of the method for first 
order differential equations with variable coefficients, that is 



(64) /' (x) + a (x) f (x) = 0, 
having the exact solution 

Xl 

(65) /(a; 2 )=exp| a (x) dx I / (xi) . 
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In this case, the only root of (jHJ) is k{x) — —ao(x), all matrices become scalar with 
D = 1 and C = 0, and the transfer exponent reduces to 



(66) M Xl -> X2 = J x-^- [do (x)} dx = x 2 a (x 2 ) - x x a (x{) - J a (x) dx. 

xi xi 

Clearly (|43[l holds and therefore, accordingly one has 



(67) F(x 2 )=eyip 
From ©, we have 



■1-2 

x 2 a (x 2 ) - xia a (xi) - J a (x) 



dx 



F( Xl ). 



(68) / (x) = exp [xk (x)} F (x) . 
Therefore, using (|67|) and i|68|) one can obtain l|65|l . 

5.2. Second-order ODEs. In this section, we first consider the simplest second- 
order differential equation with variable coefficients, given by 

(69) /" (x) + a (x) f (x) = 0. 

This equation has been studied since 1836 by Sturm [42] and has been considered 
in many reports [43-48]. Actually, any second-order ODE can be rewritten in the 
form of (|69|) by the transformation given in l|53|l . Also, the first order non- linear 
Riccati equation [2] takes the above form after suitable transformations. 

In physics this equation models the propagation of transverse-electric (TE) po- 
larized light in one-dimensional (ID) non- homogeneous media (with a position- 
dependent refractive index), or motion of a single electron in a ID potential trap. 
In optical problems / (x) is the amplitude of transverse electric field and ao (x) = 
k^tr (x) — /3 2 , in which fco is the wavenumber of the radiation in vacuum, e r (x) 
is the relative permittivity function of the medium, and (3 is the propagation 
eigenvalue. In quantum mechanics / (x) is the probability wave function and 
ao (x) = 2m [E — V {x)]/h 2 , where m is the electron mass, V (x) is the electric 
potential, H is the reduced Planck constant, and E is the energy. Radial functions 
of axial gravity waves of a black hole are governed by a generalized Klein-Gordon 
or Regge- Wheeler equation [49,50], which takes a similar form to (|69|l . 

Here, we consider the solution of this equation by DTMM. For n — 2, the kernel 
matrix U (x) can be simplified from (|36|l to 



(70) V(x) 



T-^j-exp [+x(ki - k 2 )\ 



I i_exp [-x{h-k 2 ) 



From ©, one has kf{x) + ao(x) = 0, i = 1, 2, that is k\ (x) = —jk (x) and k 2 (x) 
+jk (x), with k (x) = a/ ao (x). Therefore Ij70|l simplifies as 

k' (x) 



(71) 



XJ(x) 



2k (x) 



— l+j2k(x)x exp [j 2x k (x)] 
exp [— 2jxk (x)] —l—j2k{x)x 
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which is in agreement with our previous results [18,19,22]. The transfer matrix 
Q Xl ^ l2 is then found by (|42J) (exact evaluation of matrix exponential in this case is 
possible by the theorem described in the Appendix B). The determinant of Q Xl -*x 2 
is 



, 79 x , n , _ fc 2 (xi) - fci (xi) _ fc(xi) 

( ] n^i- k2 (x 2 )-k 1 (x 2 )- k( X2 y 

The complete solution of 169(1 is then given by J^J as 



(73) / (x) = fx (x) exp [-jxk (x)] + f 2 (x) exp [+jxk (x)] . 

The singularities of 1(691) correspond to the turning points at which both wavenum- 
bers become zero. In fact, a singularity separates the regions in which the waves 
are evanescent and propagating, or in other words, ao (x) is respectively positive 
or negative. The treatment method for singularities as described in sub-section 4.5 
can be applied to find the transfer matrix over the singularity, being approximated 
by its corresponding jump transfer matrix. 

We can classify the singularities depending on the nature of the waves across 
the singularity. Here, a singularity at x = £ is characterized as /ci(£) = k 2 (£), or 
a o(£) = 0- Correspondingly, the singularity at x — £ is referred to as type A if 
ao (x > £) > and ao (x < £) < 0, type B if ao (x > £) < and ao (x < £) > 0, and 
otherwise type C. Due to 1(16(1 the jump transfer matrix from region 1 to 2 across 
the interface x = £ for ^ takes the form [18,19,22,51] 



(74) 



Q 



l->2 



fc2 + fci r +i£(k 9 -ki) 
fc2-fci p -jg(fc2+fci ) 



2fc 2 

2/c 2 e 



2fe 2 

Here, k\ = — ao(£ — Sx) and k\ = — ao(£ + £x). Therefore, the jump transfer matrix 
over singularities given by Q^-s x -*^+Sx simplifies as 



(75) 
(76) 
(77) 



Q,£-8x-^>£+5x = 
5x— >£-\-5x 



1+1 kzl 

2 2 
2 2 



1-j 1+2 

2 2 

1+3 1-j 

2 2 



1 
1 



, typeA 



typeB 



typeC. 



The results here are obtained by setting k (£ ± Sx) fts yj±8xa' (| ± Sx) through 
the corresponding Taylor's expansion of ao (x) around x = £ and taking the limit 
-► 0+. 

The DTMM has been used to solve l(69|) numerically [18,19], and the results have 
been compared to other approaches. In cases that analytical solutions are known, 
the results have been in agreement. While the exact transfer matrix is given by 
1(42(1 . we have noticed that the reduced form as in 1(43(1 can be also used and leads 
to accurate solutions. Also, energy dispersion of an electron in the ID potential 
of infinite monatomic chain has been computed [22] by this approach. In general, 
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computation times for DTMM are at least one order of magnitude lower than the 
other existing approaches. We also have discussed how it is possible to extend 
this approach to transverse-magnetic (TM) polarized modes and inhomogeneous 
anisotropic media, where four field components are solved for at once [22]. 

We also have recently exploited the DTMM to periodic electromagnetic struc- 
tures [39], in which ao(x), and thus the wavenumber functions ki(x), i — 1, 2 are all 
periodic. This has led us to a novel yet simple mathematical description of waves 
in inhomogeneous periodic media. 

5.3. Fourth-order ODEs. As another application example, we obtain the differ- 
ential transfer matrix of a simple fourth-order differential equation that has been 
discussed in several reports [47,48] 



(78) /W (x) + a (x) f (x) = 0. 

This time from ((HJ) we have kf(x) + ao(x) = Q,i = 1, ... ,4, and thus k\ (x) = 
—k (x), k 2 (x) — —jk (x), ^3 (x) = +k (x), and k± (x) = +jk (x), in which k (x) = 
v/— do (x). The kernel matrix U (x) is found by Mathematica as 

(79) 

U(*) = |^x 

2xk(x)-3 (l+j) e (l-J>feO>0 e 2xk(x) ^ _ ^ e (l+j)xk(x) 

(1 - j) e O'-i)**0<0 2jxk (x) - 3 (1 + j) e (i+i)**(») e ^xk(x) 

e -2xk(x) ^ _ ^ e -{i+j)xk{x) _ 2xk ( x ) _ 3 (! + j) e (j-l)xk(x) 

(1 + j) e-^+J')^) e - 2jxk ^ (1 - j] \ e £-3>k(x) -2jxk(x)-3 

Here, it is possible to check the determinant of Q Xl -*x 2 given by 



/ x 2 

(80) |Q !Cl ^ X3 |=exp(tr{M :Cl ^ X3 })=expl -6 J 



k' (x) \ _ k e (gi) 

k(x) k*( X2 y 



in justification of (|52|l . 

To show the applicability of the approach, we take ao (x) — — a 4 a;~ 4 for which 
(|78|l becomes an Euler-Cauchy equation and therefore has an exact solution given 

by 



(81) / (x) = Cl x mi + c 2 x m2 + c 3 x rn3 + c 4 x mi , 

where mi = | ± ± 4-\/l + a 4 , i — 1, . . . , 4, and Ci,i — 1, . . . , 4 are arbitrary 
constants. Since k (x) = ax^ 1 , the kernel matrix (|79J) simplifies into 



(82) 
U(x) 



-l 

2* 



2a -3 (l+j)^ 1 -^ e 2a (l-j)e (1+ ^ a 

(l-j)e^- 1 ^ 2ja-3 (1 + j) e {1+ ^ a e 2ja 

(1 - j)e- {1+ ^ a -2a -3 (l+j)e { J-^ a 



(l-j)e( 1 -J') a -2aj-3 
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where N is a traceless constant matrix with eigenvalues Xi = ± ^ \/5 ± 4 + <?, 
i = 1, ... ,4. For this case, the kernel matrix U(x) commutes with the transfer 
exponent Mi^ so that (|43|1 becomes exact. By integration of (|82(l we find 

(83) Mi^ = ^I + Ina;N. 

Obviously I and N commute, so that 



(84) Q^ x = exp (Mi^ x ) = x 3/2 exp (N) = x 3/2 R [exp (A, In a;) <%] R \ 

where the constant matrix R is the diagonalizer of N. Therefore, |Qi_> x | = a; 6 , as 
required by (|80|l . Furthermore, the elements of Qi^^ involve linear combinations 
of x Ti , % = 1, . . . , 4, for which ~ rrii hold. 
Finally from ©, 0, and H2J we have 



(85) /(*)= exp [# (z)] *Q 1 ^ ar P(l) ) 

where <fr (x) has become a constant vector, and F (1) is a vector of arbitrary con- 
stants, to be determined by boundary conditions. However, comparing to H81|) 
shows that (|85|l must be indeed the general solution of i|78|) . 



5.4. Abel-Liouville-Ostogradski Formula. Found in 1827 by Abel for second- 
order differential equations and by Liouville and Ostogradsky in 1838 for the general 
case, the Wronskian should satisfy [26-28] 



(86) W 1 (x) + a n -x (x) W (x) = 0, 

which implies that the Wronskian must be essentially a constant if a„_i (x) = 0. 
Now we show that this equation can be readily reconstructed by DTMM. 

Following the discussions in sub-section 4.4 the Wronskian determinant takes 
the simple form W = |D| |exp (xK)| IQc^^l |F|, where c G S is a constant and F 
is the matrix of independent vectors as defined in section 4.4. Therefore using fTSJl 
and H51|) we find 



W = II ih ( x ) - k j ( x )] ex P i xtl {K}] 



i>j 



exp 



(87) 



in-i (c) - xa n -i (x) - J a n -i (x) dx 



exp [ca n -i (c)] J] [ki (c) - kj (c)] |F| exp 

i>3 



TT fcj(c)-fcj(c) 

11 ki(x)-kj(x) 
t>j 
x 

J a„_i (x)dx 



A exp 



/ a n _i (x) dx 



in which A is a constant. Clearly the Wronskian given by (|87(l satisfies the differ- 
ential equation JHBJ). 
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6. Conclusions 

We have presented a new analytical formalism for solution of linear homogeneous 
ordinary differential equations with variable coefficients. The formalism is based on 
the definition of jump transfer matrices and their extension into differential form. 
We presented the derivative lemma and the fundamental theorem of differential 
transfer matrix method, which support the exactness of the formalism. We also 
discussed the linear independency of solutions. The method is completely general, 
but fails in the presence of identical roots of characteristic equation. We have dis- 
cussed a method to deal with corresponding singularities. The main advantage of 
the presented method is that it deals directly with the evolution of linear indepen- 
dent solutions, rather than their derivatives. The DTMM, as described in §5, when 
applied to wave propagation problems, leads to a novel approach for understanding 
the behavior of physical systems. 

Acknowledgement. The authors wish to thank Prof. Beresford Parlett and Dr. 
Michael Parks at University of California at Berkeley for valuable discussions. 



Appendix A. Proof of (|4*£|) 

We prove the validity of l|48|) . by showing first that both sides have the same 
derivative. From [39] the left-hand-side of (|48|l should obey the differential equation 



d 
dx 2 



exp 



X2 X"2 

/ H(t) _1 H' (t)dt = exp J Hfij-'H' (t)dt 



tr {H- 1 (ar 2 ) H' (x 2 )} 



Similarly, application of the derivative theorem for determinants [32, p. 178] to the 
right-hand-side of ll4*8|) results in 



(89) 



— \U- 1 (X 1 )U(X 2 )\ = .„ 

dx2 |H (xi 



I d|H(a 2 )| |H(n 2 ) 



ti-{n-\x 2 )n'(x 2 )} , 



dx 2 |H(xi _ 

Therefore, both sides satisfy the same differential equation, and thus 



(90) 



exp 



X 2 

J H(a:) _1 H' {x)dx 



= w{x l )\R- l {x l )n(x 2 )\ 



where u>(-) is a function of only x\. But upon setting x\ = x 2 in 19U|) . one gets 
w(-) = 1 and hence the claim. 



Appendix B. Matrix Exponential of 2x2 Matrices 

The difficulties in evaluation of matrix exponential have been well known for a 
long time [52,53]. While computation of matrix exponential is possible by general 
[54-56] and approximate [57] methods, simple analytical expression exists only for 
few cases including the Euler-Rodrigues formula [50] for 3x3 skew-symmetric 
matrices. Here, we report an exact expansion of matrix exponential for arbitrary 
2x2 square matrices. 
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Assume that M is a square matrix with arbitrary complex elements. We de- 
fine the traceless matrix A = M ^tr{M}I. Since A and itr{M}I obviously 
commute, one can deduce that 

(91) exp(M) = cxp(±tr{M})cxp(A). 
Now, A satisfies the property A 2 = |A| I, so that [58] 

(92) exp (A) = cos (S) I + sine (5) A, 
in which 6 = \J— |A| and sinc(x) = sin(x)/x. Finally, 

(93) exp (M) = exp (±tr {M}) (cos 51 + sinc<5 A) . 
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